%% load data
clear; clc; close all force;
load stream.mat;

savefigures = true;
statistics = true;

alpha = 1;
corr_scatter = 150;
fontsize = 32;
lw = [1 3];
main_marker_size = 15;
markers = {'o','+'};
scatter_size = 60;

%% figure 1: stimulus for tonal stream segregation

% prepare figure
f = figure(1); clf; hold on;
f.WindowState = 'maximized';
colors = get(gca,'ColorOrder');

% tiled layout
tiledlayout(2,1,'TileSpacing','compact');
hold on;

% draw jitter bars
plot(2.8*[1 1],[0 2.2],'k--','linewidth',2)
plot(2.7*[1 1],[2.2 2.5],'k','linewidth',2)
plot(2.9*[1 1],[2.2 2.5],'k','linewidth',2)
plot([2.7 2.9],[2.55 2.55],'k:+','linewidth',2)
text(2.8,2.6,'Jitter','fontsize',24,'HorizontalAlignment','center')

% draw stimulus intervals
for a = 1:8

    % regular tempo
    start = (a-1) * 0.4;
    stop = start + 0.1;
    patch([start,start,stop,stop],[0,1,1,0],colors(1,:),'FaceAlpha',1/2);

    % jittered
    start = start + 0.2 + 0.1*(rand-0.5);
    stop = start + 0.1;
    patch([start,start,stop,stop],[0,1,1,0]+1,colors(2,:));

end

% axes
set(gca,'fontsize',24,'xlim',[0 3.2],'xtick',0:0.4:2.8,'ylim',[0 2.7],'ytick',[0.5 1.5],'yticklabel',{'Target','Masker'})
xlabel('Time (seconds)')

% save figure
if savefigures
    drawnow; pause(1/4);
    saveas(f,'figures/fig 1 stream stimulus.png')
end

%% figure 2: stimulus for speech segregation 

% prepare figure
f = figure(2); clf; hold on;  
f.WindowState = 'maximized';
colors = get(gca,'ColorOrder');

% stop time
stop = 2;

% tiled layout
tiledlayout(2,1,'TileSpacing','compact');

% sparse noise temporal envelope
nexttile(1); hold all;

% sparse noise
[sparse,fs] = audioread('crisp_babble_sparse_pitchshift_0.wav');
time = (1:size(sparse,1)) / fs;
plot(time(time<stop),sparse(time<stop))

% target
[x,fs] = audioread('AIRPLANE.wav');
time = (1:size(x,1)) / fs;
plot(1+time(time<stop),x(time<stop),'color',[colors(2,:),1/2])

% axis
set(allchild(gca),'linewidth',2)
set(gca,'FontSize',24,'linewidth',2,'xticklabel','','ylim',0.8*[-1 1])
title('Sparse Noise')
ylabel('Amplitude')

% dense noise temporal envelope
nexttile(2); hold all;

% dense noise
[dense,fs] = audioread('crisp_babble.wav');
time = (1:size(dense,1)) / fs;
plot(time(time<stop),dense(time<stop))

% target
[x,fs] = audioread('AIRPLANE.wav');
time = (1:size(x,1)) / fs;
plot(1+time(time<stop),x(time<stop),'color',[colors(2,:),1/2])

% axes
set(allchild(gca),'linewidth',2)
set(gca,'FontSize',24,'linewidth',2,'ylim',0.8*[-1 1])
title('Dense Noise')
xlabel('Time (s)')
ylabel('Amplitude')

% % compare power spectral densities
% nexttile(2,[2,1]); hold all;
% 
% % power spectral density (sparse)
% S = db(abs(fft(sparse)));
% f = (1:size(S,1)) * (fs / size(S,1));
% plot(f(f<10e3),S(f<10e3)+30)
% 
% % power spectral density (dense)
% D = db(abs(fft(dense)));
% f = (1:size(D,1)) * (fs / size(D,1));
% plot(f(f<10e3),D(f<10e3)+20)
% 
% % axes
% set(gca,'FontSize',24,'ylim',[0 100])

% save figure
if savefigures
    drawnow; pause(1/4);
    saveas(f,'figures/fig 2 speech stimulus.png')
end

%% figure 3: frequency and f0 discrimination

% figure
f = figure(3); clf; hold on;
f.WindowState = 'maximized';
colors = get(gca,'ColorOrder');

% tiled layout
tiledlayout(1,3,'TileSpacing','compact');

% Pure Tones
nexttile(1); hold on;

% loop through groups
for group = 1:2

    % switch group
    switch group
        case 1
            data = nanmean(NH_FDT,2);
            I = data > 1.5; data(I) = [];
        case 2
            data = nanmean(CI_FDT,2);
    end

    % abscissa positioning
    x = 1+(group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,120,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/4)

    % bootstrapped confidence intervals
    m = bootstrp(1000,@mean,data);
    [fi,xi] = ksdensity(m);
    c = cumsum(fi)/sum(fi);

    % 99% Confidence Intervals
    interval = c > 0.01 & c < 0.99;
    X1 = x + fi(interval)/max(fi(interval))/8;
    X2 = x - fi(interval)/max(fi(interval))/8;
    Y1 = xi(interval);
    plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)

    % 95% Confidence Interval
    interval = c > 0.05 & c < 0.95;
    X1 = x + fi(interval)/max(fi(interval))/8;
    X2 = x - fi(interval)/max(fi(interval))/8;
    Y1 = xi(interval);
    L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    % means and standard errors
    errorbar(x,mean(data),std(data,[],1)/sqrt(size(data,1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',12)

end

% axes
set(findobj(gca,'type','line'),'linew',2)
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 1.5],'xtick',1,'xticklabel',{'1000'},...
    'ylim',[-1.2 log10(300)],'ytick',-1:2,'yticklabel',{'0.1','1','10','100'});
title('Pure Tones')
xlabel('Frequency (Hz)')
ylabel({'Discrimination Threshold';'(% difference)'})

% right axes
yyaxis right
set(gca,'fontsize',fontsize,...
    'ylim',[-1.2 log10(300)],'ytick',log10(100*(2.^([1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16]/12)-1)),...
    'yticklabel',{'','','','','','','','','','',''});
ax = gca;
ax.YAxis(1).Color = 'k';
ax.YAxis(2).Color = 'k';

% 3-semitone line
plot([0.5 1.5],[1 1]*log10(100*(2^(3/12)-1)),'k--','LineWidth',2);

% Complex Tones
nexttile(2); hold on

% loop through groups
for group = 1:2
    
    % switch group
    switch group
        case 1
            data = nanmean(NH_F0DT,2);
            I = data > 1.5; data(I) = [];
        case 2
            data = nanmean(CI_F0DT,2);
    end

    % abscissa position
    x = 1+(group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,120,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/4)

    % bootstrapped confidence intervals
    m = bootstrp(1000,@mean,data);
    [fi,xi] = ksdensity(m);
    c = cumsum(fi)/sum(fi);

    % 99% Confidence Intervals
    interval = c > 0.01 & c < 0.99;
    X1 = x + fi(interval)/max(fi(interval))/10;
    X2 = x - fi(interval)/max(fi(interval))/10;
    Y1 = xi(interval);
    plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)

    % 95% Confidence Interval
    interval = c > 0.05 & c < 0.95;
    X1 = x + fi(interval)/max(fi(interval))/10;
    X2 = x - fi(interval)/max(fi(interval))/10;
    Y1 = xi(interval);
    L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    % means and standard erros
    errorbar(x,mean(data),std(data,[],1)/sqrt(size(data,1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',12)

end

% axes
set(findobj(gca,'type','line'),'linew',2)
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 1.5],'xtick',1:2,'xticklabel',{'110'},...
    'ylim',[-1.2 log10(300)],'ytick',-1:2,'yticklabel',{'','','',''});
title('Complex Tones')

% right axes
yyaxis right
set(gca,'fontsize',fontsize,...
'ylim',[-1.2 log10(300)],'ytick',log10(100*(2.^([1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16]/12)-1)),...
'yticklabel',{'1/32','','1/16','','1/4','','1','','4','','16'});
xlabel('Fundamental Frequency (Hz)')
ylabel({'(semitones)'})
ax = gca;
ax.YAxis(1).Color = 'k';
ax.YAxis(2).Color = 'k';

% 3-semitone line
plot([0.5 1.5],[1 1]*log10(100*(2^(3/12)-1)),'k--','LineWidth',2);

% legend
h(1) = plot(0,0,'o','Color',colors(1,:),'markersize',14,'MarkerFaceColor',colors(1,:));
h(2) = plot(0,0,'+','Color',colors(2,:),'linewidth',3,'markersize',14); 
legend(h,{'No known hearing loss','Cochlear implant users'},'Box','off','Location','southwest','fontsize',20);

% statistics
if statistics

    % NH_FDT is subject x repetition
    % NH_F0DT is subject x repetition

    %
    disp(' --- FDTs --- ')

    % means
    x = 10^mean(mean(NH_FDT)); y = 10^mean(mean(CI_FDT)); 
    disp(['NH FDT: ' num2str(x,2) '%, ' num2str(12*log2(1+x/100),2) ' semitones'])
    disp(['CI FDT: ' num2str(y,2) '%, ' num2str(12*log2(1+y/100),2) ' semitones'])

    % compare NH and CI FDTs
    x = mean(NH_FDT,2); 
    y = mean(CI_FDT,2);
    m1 = mean(x);
    m2 = mean(y);
    v1 = var(x);
    v2 = var(y);
    cohen = (m2-m1)/sqrt(v1+v2); disp(['Cohen: ' num2str(cohen)]);
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);

    % F0DTs
    disp(' --- F0DTs --- ')

    % means
    x = 10^mean(mean(NH_F0DT)); y = 10^mean(mean(CI_F0DT)); 
    disp(['NH F0DT: ' num2str(x,2) '%, ' num2str(12*log2(1+x/100),2) ' semitones'])
    disp(['CI F0DT: ' num2str(y,2) '%, ' num2str(12*log2(1+y/100),2) ' semitones'])

    % compare NH and CI F0DTs
    x = mean(NH_F0DT,2); 
    y = mean(CI_F0DT,2);
    m1 = mean(x);
    m2 = mean(y);
    v1 = var(x);
    v2 = var(y);
    cohen = (m2-m1)/sqrt(v1+v2); disp(['Cohen: ' num2str(cohen)]);
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);

    % compare FDTs to F0DTs
    disp(' --- FDTs vs F0DTs --- ')

    % 
    disp(' --- NH --- ')
    x = mean(NH_F0DT,2); 
    y = mean(NH_FDT,2);
    m1 = mean(x);
    m2 = mean(y);
    v1 = var(x);
    v2 = var(y);
    cohen = (m2-m1)/sqrt(v1+v2); disp(['Cohen: ' num2str(cohen)]);
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);

    % 
    disp(' --- CI --- ')
    x = mean(CI_F0DT,2); 
    y = mean(CI_FDT,2);
    m1 = mean(x);
    m2 = mean(y);
    v1 = var(x);
    v2 = var(y);
    cohen = (m2-m1)/sqrt(v1+v2); disp(['Cohen: ' num2str(cohen)]);
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);

    % organize data for ANOVA
    data = [
        NH_FDT, NH_F0DT;
        CI_FDT, CI_F0DT;
    ];

    % between factors
    group = [ones(size(NH_FDT,1),1);2*ones(size(CI_FDT,1),1)];

     % measurement table
    t = array2table([data group]);
    t.Properties.VariableNames(end) = {'group'};
    t.group = categorical(t.group);

    % within factors
    stimulus = []; repetition = [];
    for a = 1:2 % stimulus
        for b = 1:3 % repetition
            stimulus(end+1) = a;
            repetition(end+1) = b;
        end
    end

    % within design
    within = table(stimulus',repetition','VariableNames',{'stimulus','repetition'});
    within.stimulus = categorical(within.stimulus);
    within.repetition = categorical(within.repetition);

    % fit the model
    rm = fitrm(t,'Var1-Var6~group','WithinDesign',within);

    % run the analyses
    anova_resolution = ranova(rm,'WithinModel','stimulus');
    mc_resolution = multcompare(rm,'stimulus','By','group','ComparisonType','bonferroni');

    % add statistics to plot
    A = anova_resolution;
    nexttile(3); hold on; axis off; fs = 22; y_pos = 1; o1 = 0.05; o2 = 0.08; ind = 2;
    text(0,y_pos,'\bf{group} ','FontSize',fs); y_pos = y_pos - o1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group x stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' sprintf('%.2f',A.pValue(ind))],'FontSize',fs);
  
end

% save figure
if savefigures
    drawnow; pause(1/4);
    saveas(f,'figures/fig 3 discrimination.png')
end

%% figure 4: tonal stream segregation

% prepare figure
f = figure(4); clf; hold on;
f.WindowState = 'maximized';
colors = get(gca,'ColorOrder');
lw = [1,2];

% tiled layout
tiledlayout(1,5,'TileSpacing','compact');
nexttile(1,[1 2]);hold on

% loop through groups
for group = 1:2

    % switch group
    switch group
        case 1
            data = mean(NH_SS_F,3);
        case 2
            data = mean(CI_SS_F,3);
    end

    % abscissa positioning
    x = (1:5)+(group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,60,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/8);

    % bootstrapped confidence intervals
    for condition = 1:5
        
        % bootstrapped confidence intervals
        m = bootstrp(1000,@mean,data(:,condition));
        [fi,xi] = ksdensity(m);
        c = cumsum(fi)/sum(fi);
    
        % 99% Confidence Intervals
        interval = c > 0.01 & c < 0.99;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)
    
        % 95% Confidence Interval
        interval = c > 0.05 & c < 0.95;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    end

    % means and standard erros
    errorbar(x,mean(data),std(data,[],1)/sqrt(size(data,1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',8)

end

% axes
ytick = [1:10,20:10:100];
yticklabel = strings(size(ytick));
yticklabel(1) = '1';
yticklabel(10) = '10';
yticklabel(19) = '100';
set(findobj(gca,'type','line'),'linew',2)
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 5.5],'xtick',1:5,'xticklabel',{'0','3','6','12','Quiet'},'XTickLabelRotation',0,...
    'ylim',[log10(3) log10(200)],'ytick',log10(ytick),'yticklabel',yticklabel);
title('Pure Tones')
xlabel('Frequency Distance (semitones)');
ylabel({'Jitter Detection Threshold (ms)'})

% ceiling
plot([.5 5.5],[2 2],'k-','linewidth',2);
text(.6,2.12,'ceiling','fontsize',30);

% Stream Segregation (Complex Tones)
nexttile(3,[1 2]); hold on;

% loop through groups
L  = [];
for group = 1:2
   
    % switch group
    switch group
        case 1
            data = mean(NH_SS_F0,3);
        case 2
            data = mean(CI_SS_F0,3);
    end

    % abscissa positioning
    x = (1:5)+(group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,60,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/8);

    % bootstrapped confidence intervals
    for condition = 1:5

        % bootstrapped confidence intervals
        d = data(:,condition);
        m = bootstrp(1000,@mean,d(d>0));
        [fi,xi] = ksdensity(m);
        c = cumsum(fi)/sum(fi);
    
        % 99% Confidence Intervals
        interval = c > 0.01 & c < 0.99;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)
    
        % 95% Confidence Interval
        interval = c > 0.05 & c < 0.95;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    end

    % sample means and standard errors (note, data > 0 correction for one sample)
    data(data==-Inf) = nan;
    errorbar(x,nanmean(data),std(data(~isnan(data)),[],1)/sqrt(size(data(data>0),1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',8)

end

% axes
set(findobj(gca,'type','line'),'linew',2)
ytick = [1:10,20:10:100]; 
yticklabel = strings(size(ytick));
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 5.5],'xtick',1:5,'xticklabel',{'0','3','6','12','Quiet'},'XTickLabelRotation',0,...
    'ylim',[log10(3) log10(200)],'ytick',log10(ytick),'yticklabel',yticklabel);
title('Complex Tones');
xlabel('F0 Distance (semitones)');

% ceiling
plot([.5 5.5],[2 2],'k-','linewidth',2);
text(.6,2.12,'ceiling','fontsize',30);

% legend
h(1) = plot(0,0,'o','Color',colors(1,:),'markersize',14,'MarkerFaceColor',colors(1,:)); 
h(2) = plot(0,0,'+','Color',colors(2,:),'linewidth',3,'markersize',14); 
legend(h,{'No known hearing loss','Cochlear implant users'},'Box','off','Location','southwest','fontsize',20);

% cohens d
x = mean(NH_SS_F0(:,5,:),3); 
y = mean(NH_SS_F0(:,4,:),3);
m1 = mean(x);
m2 = mean(y);
v1 = var(x);
v2 = var(y);
cohen = (m2-m1)/sqrt(v1+v2); disp(cohen);

% statistics
if statistics

    % NH_SS_F is subject x snr x repetition
    % NH_SS_F0 is subject x snr x repetition

    % organize data
    data = [
        reshape(NH_SS_F,size(NH_SS_F,1),numel(NH_SS_F)/size(NH_SS_F,1)),reshape(NH_SS_F0,size(NH_SS_F0,1),numel(NH_SS_F0)/size(NH_SS_F0,1));
        reshape(CI_SS_F,size(CI_SS_F,1),numel(CI_SS_F)/size(CI_SS_F,1)),reshape(CI_SS_F0,size(CI_SS_F0,1),numel(CI_SS_F0)/size(CI_SS_F0,1))
        ];

    % replace one abberant value
    data(data<0) = 2;

    % grand average cohen's d
    m1 = mean(mean(data(group==1,:),2));
    v1 = var(mean(data(group==1,:),2));
    m2 = mean(mean(data(group==2,:),2));
    v2 = var(mean(data(group==2,:),2));
    cohen1 = (m2 - m1) / sqrt(v1+v2); disp(cohen);

    % between factors
    group = [ones(size(NH_SS_F,1),1);2*ones(size(CI_SS_F,1),1)];

    % measurement table
    t = array2table([data group]);
    t.Properties.VariableNames(end) = {'group'};
    t.group = categorical(t.group);

    % within factors
    stimulus = []; snr = []; repetition = [];
    for a = 1:2 % stimulus
        for b = 1:5 % snr
            for c = 1:3 % repetition
                stimulus(end+1) = a;
                snr(end+1) = b;
                repetition(end+1) = c;
            end
        end
    end

    % within design
    within = table(stimulus',snr',repetition','VariableNames',{'stimulus','pitch','repetition'});
    within.stimulus = categorical(within.stimulus);
    within.pitch = categorical(within.pitch);
    within.repetition = categorical(within.repetition);

    % fit the model
    rm = fitrm(t,'Var1-Var30~group','WithinDesign',within);

    % run the analyses
    anova_stream = ranova(rm,'WithinModel','stimulus*pitch');
    mc_stream = multcompare(rm,'stimulus','By','group','ComparisonType','bonferroni');
    
    % add statistics to plot
    A = anova_stream;
    nexttile(5,[1 1]); hold on; axis off; fs = 22; y_pos = 1; o1 = 0.05; o2 = 0.08; ind = 2;
    text(0,y_pos,'\bf{group} ','FontSize',fs); y_pos = y_pos - o1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group x stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{frequency distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group x frequency distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{stimulus x frequency distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' sprintf('%.2f',A.pValue(ind))],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,'\bf{x stimulus x frequency distance} ','FontSize',fs); y_pos = y_pos - o1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' sprintf('%.2f',A.pValue(ind))],'FontSize',fs);
    
end

% save figure
if savefigures
    drawnow; pause(1);
    saveas(f,'figures/fig 4 stream segregation.png')
end
%% figure 5: speech reception thresholds

% figure
f = figure(5); clf; hold on;
f.WindowState = 'maximized';
colors = get(gca,'ColorOrder');
lw = [1,2];

% sparse noise
tiledlayout(1,5,'TileSpacing','compact');
nexttile(1,[1 2]); hold on;

% loop through groups
for group = 1:2

    % specify data based on group
    switch group
        case 1
            data = mean(NH_SR_Sparse,3);
        case 2
            data = mean(CI_SR_Sparse,3);
    end

    % Cohen's d
    x = data(:,1); y = data(:,2);
    m1 = mean(x); m2 = mean(y);
    v1 = var(x); v2 = var(y);
    d = (m2 - m1) / sqrt(v1+v2);
    disp(['Sparse Noise, Group ' num2str(group)])
    disp([num2str(m1,3),', ',num2str(m2,3),', ',num2str(m1-m2,3),', ',num2str(d,3)])
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);
    disp(' ');

    % abscissa
    x = (1:5) + (group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,60,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/4);

    % bootstrapped confidence intervals
    for condition = 1:5

        % bootstrapped confidence intervals
        d = data(:,condition);
        m = bootstrp(1000,@mean,d);
        [fi,xi] = ksdensity(m);
        c = cumsum(fi)/sum(fi);
    
        % 99% Confidence Intervals
        interval = c > 0.01 & c < 0.99;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)
    
        % 95% Confidence Interval
        interval = c > 0.05 & c < 0.95;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    end

    % means and standard erros
    errorbar(x,nanmean(data),std(data(~isnan(data)),[],1)/sqrt(size(data,1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',8)

end

% axes
set(findobj(gca,'type','line'),'linew',2)
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 5.5],'xtick',1:5,'xticklabel',{'0','3','6','12','Quiet'},'XTickLabelRotation',0,...
    'ylim',[-60 20],'ytick',-60:20:20,'yticklabel',{'-60','-40','-20','0','20'});
title('Sparse Noise')
ylabel('Speech Reception Threshold (dB SNR)')

% Dense
nexttile(3,[1 2]); hold on

% loop through groups (plotting)
for group = 1:2

    % switch mode
    switch group
        case 1
            data = mean(NH_SR_Dense,3);
        case 2
            data = mean(CI_SR_Dense,3);
    end

    % Cohen's d
    x = data(:,1); y = data(:,2);
    m1 = mean(x); m2 = mean(y);
    v1 = var(x); v2 = var(y);
    d = (m2 - m1) / sqrt(v1+v2);
    disp(['Dense Noise, Group ' num2str(group)])
    disp([num2str(m1,3),', ',num2str(m2,3),', ',num2str(m1-m2,3),', ',num2str(d,3)])
    [~,p,~,s] = ttest2(x,y); % associated t test
    disp([ 't(' num2str(s.df) ') = ' num2str(s.tstat) ', p = ' num2str(p) ]);
    disp(' ');

    % abscissa positioning
    x = (1:5)+(group-3/2)/4;

    % individual data
    swarmchart(ones(size(data,1),1)*x,data,60,colors(group,:),'linewidth',lw(group),...
        'marker',markers{group},'markerfacecolor',colors(group,:),'xjitterwidth',1/4);

    % bootstrapped confidence intervals
    for condition = 1:5

         % bootstrapped confidence intervals
        d = data(:,condition);
        m = bootstrp(1000,@mean,d(~isnan(d)));
        [fi,xi] = ksdensity(m);
        c = cumsum(fi)/sum(fi);
    
        % 99% Confidence Intervals
        interval = c > 0.01 & c < 0.99;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        plot([X1 fliplr(X2) X1(1)],[Y1 fliplr(Y1) Y1(1)],'color',colors(group,:),'linewidth',2)
    
        % 95% Confidence Interval
        interval = c > 0.05 & c < 0.95;
        X1 = x(condition) + fi(interval)/max(fi(interval))/8;
        X2 = x(condition) - fi(interval)/max(fi(interval))/8;
        Y1 = xi(interval);
        L(group) = patch([X1,fliplr(X2)],[Y1,fliplr(Y1)],colors(group,:),'FaceAlpha',1/2);

    end

    % means and standard erros
    errorbar(x,nanmean(data),std(data(~isnan(data)),[],1)/sqrt(size(data,1)),'color','k','linestyle','none','linewidth',2,...
        'marker',markers{group},'markerfacecolor','k','markersize',8)

end

% axes
set(findobj(gca,'type','line'),'linew',2)
set(gca,'fontsize',fontsize,'box','off','linewidth',2,...
    'xlim',[0.5 5.5],'xtick',1:5,'xticklabel',{'0','3','6','12','Quiet'},'XTickLabelRotation',0,...
    'ylim',[-60 20],'ytick',-60:20:20,'yticklabel',{'','','','',''});
title('Dense Noise')

%  legend
legend(L,{'No known hearing loss','Cochlear implant users'},'fontsize',20,'Location','southwest');
legend box off;

% xlabel
h = xlabel('F0 Distance (semitones)');
set(h,'position',get(h,'position')-[2.8 -2 0])

% legend
h(1) = plot(0,0,'o','Color',colors(1,:),'markersize',14,'MarkerFaceColor',colors(1,:)); 
h(2) = plot(0,0,'+','Color',colors(2,:),'linewidth',3,'markersize',14); 
legend(h,{'No known hearing loss','Cochlear implant users'},'Box','off','Location','southwest','fontsize',20);

% statistics
nexttile(5); hold on

% statistics
if statistics

    % data
    NH_spin = [NH_SR_Sparse NH_SR_Dense];
    NH_spin = NH_spin(:,[2 3 4 5 7 8 9 10]);
    CI_spin = [CI_SR_Sparse CI_SR_Dense];
    CI_spin = CI_spin(:,[2 3 4 5 7 8 9 10]);
    data = [NH_spin;CI_spin];
    
    % group index
    group = [ones(size(NH_spin,1),1);2*ones(size(CI_spin,1),1)];
    
    % measurement index
    sparse_or_dense = [ones(1,4) 2*ones(1,4)];
    
    % condition index
    pitch_distance = repmat([1 2 3 4],1,2);
    
    % measurement table
    t = array2table([data group]);
    t.Properties.VariableNames(end) = {'group'};
    t.group = categorical(t.group);
    
    % within design
    within = table(sparse_or_dense',pitch_distance','VariableNames',{'stimulus','pitch'});
    within.stimulus = categorical(within.stimulus);
    within.pitch = categorical(within.pitch);
    
    % fit the model
    rm = fitrm(t,'Var1-Var8~group','WithinDesign',within);
    
    % run the analyses
    anova_speech = ranova(rm,'WithinModel','stimulus*pitch');
    mc_speech = multcompare(rm,'stimulus','By','group','ComparisonType','bonferroni');

    % add statistics to plot
    A = anova_speech;
    nexttile(5,[1 1]); hold on; axis off; fs = 22; y_pos = 1; o1 = 0.05; o2 = 0.08; ind = 2;
    text(0,y_pos,'\bf{group} ','FontSize',fs); y_pos = y_pos - o1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group x stimulus} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{F0 distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group x F0 distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{stimulus x F0 distance} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 2;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+2)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',22); y_pos = y_pos - o2;
    text(0,y_pos,'\bf{group} ','FontSize',fs); y_pos = y_pos - o1; ind = ind + 1;
    text(0,y_pos,'\bf{x stimulus x F0 distance} ','FontSize',fs); y_pos = y_pos - o1;
    text(0,y_pos,['{\it F}(' num2str(A.DF(ind)) ',' num2str(A.DF(ind+1)) ') = ' num2str(A.F(ind),2) ', {\it p} = ' num2str(A.pValue(ind),2)],'FontSize',fs); y_pos = y_pos - o2;
    
end

% save figure
if savefigures
    drawnow; pause(1);
    saveas(f,'figures/fig 5 speech reception.png')
end

%% figure 6: correlations between tonal streaming and discrimination
% CI_FDT: subject x repetition
% CI_F0DT: subject x repetition
% CI_SS_F: subject x condition x repetition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SS_F0: subject x condition x repetition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SR_Sparse: subject x condition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SR_Dense: subject x condition ( condition order: 0, 3, 6, 12, Quiet )

% figure
f = figure(6); clf; hold all;
f.WindowState = 'maximized';
colors = get(groot,'defaultAxesColorOrder');
groupnames = {'No Known Hearing Loss','Cochlear Implant'};

% tiled layout
tiledlayout(1,2,'TileSpacing','compact')

%
for comparison = 1:2

    %
    nexttile; hold all;

    % loop over participant groups
    for group = 1:2
    
        %
        ind = 2:4;
        switch group
            case 1

                switch comparison
                    case 1

                        x = mean(NH_FDT,2);
                        y = mean(mean(NH_SS_F(:,ind,:),3),2);

                    case 2

                        x = mean(NH_F0DT,2);
                        y = mean(mean(NH_SS_F0(:,ind,:),3),2);

                end

                % remove outlier
                 I = x > 1.5; x(I) = []; y(I) = [];
    
                %
                xt = -1; yt = 0.7; offset = 1/16;
    
            case 2

               switch comparison
                    case 1

                        x = mean(CI_FDT,2);
                        y = mean(mean(CI_SS_F(:,ind,:),3),2);

                    case 2

                        x = mean(CI_F0DT,2);
                        y = mean(mean(CI_SS_F0(:,ind,:),3),2);

                end

                    % remove outlier
                     I = y < 0; x(I) = []; y(I) = [];
    
                    %
                    xt = 1.5; yt = 1.5; offset = 1/16;
    
        end
    
        % individual data
        if group == 1, lw = 1; else, lw = 3; end
        plot(x,y,'color',colors(group,:),'linestyle','none','linewidth',lw,...
            'marker',markers{group},'markerfacecolor',colors(group,:),'markersize',12)
        
        % linear regression
        p = polyfit(x,y,1);
        x0 = [min(x) max(x)];
        y0 = polyval(p,x0);
        plot(x0,y0,'color',colors(group,:),'LineWidth',4)
        
        % correlation statistics
        [r,p] = corrcoef(x,y);
        text(xt,yt,groupnames{group},'fontsize',24,'Color',colors(group,:))
        if p(2) < 0.001
            text(xt,yt-2*offset,['{\it r} = ' num2str(r(2),2) ' (' '{\it p} < 0.001)'],'fontsize',24,'Color',colors(group,:))
        else
            text(xt,yt-offset,['{\it r} = ' num2str(r(2),2) ' (' '{\it p} = ' num2str(p(2),2) ')'],'fontsize',24,'Color',colors(group,:))
        end
        disp(num2str(r(2)^2,2))
    
    end

    % axes
    switch comparison
        case 1
            xlabel({'Frequency Discrimination','(% Difference)'})
            ylabel('Jitter Detection Threshold (ms)')
            title('Pure Tones')
        case 2
            xlabel({'F0 Discrimination','(% Difference)'})
            title('Complex Tones')
    end

    %
    ytick = [1:10,20:10:100];
    yticklabel = strings(size(ytick));
    yticklabel(1) = '1'; yticklabel(10) = '10'; yticklabel(19) = '100';
    set(gca,'fontsize',30,...
        'xlim',[-1.2 log10(300)],'xtick',-1:2,'xticklabel',{'0.1','1','10','100'},...
        'ylim',[log10(3) log10(200)],'ytick',log10(ytick),'yticklabel',yticklabel)

    % 3-semitone line
    plot([1 1]*log10(100*(2^(3/12)-1)),[log10(3) log10(200)],'k--','LineWidth',2);

end

% hack to prevent crashing on mac
drawnow; pause(1);

% save figure
if savefigures
    saveas(f,'figures/fig 6 stream correlations.png')
end
%% figure 7: correlations between speech reception and pitch resolution

% figure
f = figure(7); clf; tiledlayout(2,2);
f.WindowState = 'maximized';

%
groupnames = {'NKHL','CI Users'};

%
for comparison = 1:4

    %
    nexttile; hold all;

    %
    for group = 1:2
        disp(['group: ' num2str(group), ' comparison: ' num2str(comparison)])
    
        %
        ind = 2:4;
        switch group
            case 1
                
                switch comparison
                    case 1

                        x = mean(NH_FDT,2);
                        y = mean(nanmean(NH_SR_Sparse(:,ind,:),3),2);

                    case 2

                        x = mean(NH_F0DT,2);
                        y = mean(nanmean(NH_SR_Sparse(:,ind,:),3),2);

                    case 3

                        x = mean(NH_FDT,2);
                        y = mean(nanmean(NH_SR_Dense(:,ind,:),3),2);

                    case 4

                        x = mean(NH_F0DT,2);
                        y = mean(nanmean(NH_SR_Dense(:,ind,:),3),2);

                end

                % remove outlier
                I = x > 1.5; x(I) = []; y(I) = [];

                %
                xt = 2.1; yt = -5; offset = 5;
    
            case 2

                switch comparison
                    case 1

                        x = mean(CI_FDT,2);
                        y = mean(nanmean(CI_SR_Sparse(:,ind,:),3),2);

                    case 2

                        x = mean(CI_F0DT,2);
                        y = mean(nanmean(CI_SR_Sparse(:,ind,:),3),2);

                    case 3

                        x = mean(CI_FDT,2);
                        y = mean(nanmean(CI_SR_Dense(:,ind,:),3),2);

                    case 4

                        x = mean(CI_F0DT,2);
                        y = mean(nanmean(CI_SR_Dense(:,ind,:),3),2);

                end

                %
                xt = 2.1; yt = 15; offset = 5;
    
        end
    
        % individual data
        if group == 1, lw = 1; else, lw = 3; end
        plot(x,y,'color',colors(group,:),'linestyle','none','linewidth',lw,...
            'marker',markers{group},'markerfacecolor',colors(group,:),'markersize',12);
        
        % linear regression
        p = polyfit(x,y,1);
        x0 = [min(x) max(x)];
        y0 = polyval(p,x0);
        plot(x0,y0,'color',colors(group,:),'LineWidth',4)
        
        % correlation statistics
        [r,p] = corrcoef(x,y);
        text(xt,yt,groupnames{group},'fontsize',22,'Color',colors(group,:))
        if p(2) < 0.05 / 6
            text(xt,yt-offset,['{\bf {\it r} = ' num2str(r(2),2) '}'],'fontsize',24,'Color',colors(group,:))
            if p(2) < 0.001
                text(xt,yt-3*offset,'{\bf {\it p} < 0.001}','fontsize',24,'Color',colors(group,:))
            else
                text(xt,yt-3*offset,['{\bf {\it p} = ' num2str(p(2),2) '}'],'fontsize',24,'Color',colors(group,:))
            end
        else
            text(xt,yt-offset,['{\it r} = ' num2str(r(2),2)],'fontsize',24,'Color',colors(group,:))
            text(xt,yt-2*offset,['{\it p} = ' num2str(p(2),2)],'fontsize',24,'Color',colors(group,:))
        end
        disp(r(2)^2)
    
    end

    %
    ytick = [1:10,20:10:100];
    yticklabel = strings(size(ytick));
    yticklabel(1) = '1'; yticklabel(10) = '10'; yticklabel(19) = '100';
    set(gca,'fontsize',30,...
        'xlim',[-1.2 2.4],'xtick',-1:2,'xticklabel',{'0.1','1','10','100'},...
        'ylim',[-30 20],'ytick',-60:20:20,'yticklabel',{'','','','',''})

    %
    switch comparison
        case 1
            text(-1,20,'\bf Sparse Noise','fontsize',24)
            set(gca,'ytick',-40:20:20,'yticklabel',{'-40','-20','0','20'})
            ylabel({'Speech Reception','(dB SNR)'})
        case 2
            text(-1,20,'\bf Sparse Noise','fontsize',24)
        case 3
            text(-1,20,'\bf Dense Noise','fontsize',24)
            set(gca,'ytick',-40:20:20,'yticklabel',{'-40','-20','0','20'})
            xlabel({'Frequency Discrimination','(% Difference)'})
            ylabel({'Speech Reception','(dB SNR)'})
        case 4
            text(-1,20,'\bf Dense Noise','fontsize',24)
            xlabel({'F0 Discrimination','(% Difference)'})
    end

end

% hack to prevent crashing on mac
drawnow; pause(1);

% save figure
if savefigures
    saveas(f,'figures/fig 7 speech correlations.png')
end

%% figure 8: correlations between masking release and pitch resolution

% prepare figure
f = figure(8); clf; tiledlayout(2,2);
f.WindowState = 'maximized';

%
for comparison = 1:4
    nexttile; hold all;

    %
    for group = 1:2
    
        %
        ind = 2:4;
        switch group
            case 1
                
                switch comparison
                    case 1

                        x = mean(NH_FDT,2);
                        y = mean(nanmean(NH_SR_Sparse(:,ind,:),3),2) - nanmean(NH_SR_Sparse(:,1,:),3);
                        y = -y;

                    case 2

                        x = mean(NH_F0DT,2);
                        y = mean(nanmean(NH_SR_Sparse(:,ind,:),3),2) - nanmean(NH_SR_Sparse(:,1,:),3);
                        y = -y;

                    case 3

                        x = mean(NH_FDT,2);
                        y = mean(nanmean(NH_SR_Dense(:,ind,:),3),2) - nanmean(NH_SR_Dense(:,1,:),3);
                        y = -y;

                    case 4

                        x = mean(NH_F0DT,2);
                        y = mean(nanmean(NH_SR_Dense(:,ind,:),3),2) - nanmean(NH_SR_Dense(:,1,:),3);
                        y = -y;
                end

                % remove outlier
                I = x > 1.5; x(I) = []; y(I) = [];

                %
                xt = 1/2; yt = 30; offset = 5;
    
            case 2

                switch comparison
                    case 1

                        x = mean(CI_FDT,2);
                        y = mean(nanmean( CI_SR_Sparse(:,ind,:),3),2) - nanmean( CI_SR_Sparse(:,1,:),3);
                        y = -y;

                    case 2

                        x = mean(CI_F0DT,2);
                        y = mean(nanmean( CI_SR_Sparse(:,ind,:),3),2) - nanmean( CI_SR_Sparse(:,1,:),3);
                        y = -y;

                    case 3

                        x = mean(CI_FDT,2);
                        y = mean(nanmean( CI_SR_Dense(:,ind,:),3),2) - nanmean( CI_SR_Dense(:,1,:),3);
                        y = -y;

                    case 4

                        x = mean(CI_F0DT,2);
                        y = mean(nanmean( CI_SR_Dense(:,ind,:),3),2) - nanmean( CI_SR_Dense(:,1,:),3);
                        y = -y;

                end

                %
                xt = 1.4; yt = 20; offset = 5;
    
        end
    
        % individual data
        if group == 1, lw = 1; else, lw = 3; end
        plot(x,y,'color',colors(group,:),'linestyle','none','linewidth',lw,...
            'marker',markers{group},'markerfacecolor',colors(group,:),'markersize',12)
        
        % linear regression
        p = polyfit(x,y,1);
        x0 = [min(x) max(x)];
        y0 = polyval(p,x0);
        plot(x0,y0,'color',colors(group,:),'LineWidth',4)
        
        % correlation statistics
        [r,p] = corrcoef(x,y);
        text(xt,yt,groupnames{group},'fontsize',22,'Color',colors(group,:))
        if p(2) < 0.05 / 6
            text(xt,yt-offset,['{\bf {\it r} = ' num2str(r(2),2) '}'],'fontsize',24,'Color',colors(group,:))
            if p(2) < 0.001
                text(xt,yt-2*offset,'{\bf {\it p} < 0.001}','fontsize',24,'Color',colors(group,:))
            else
                text(xt,yt-2*offset,['{\bf {\it p} = ' num2str(p(2),2) '}'],'fontsize',24,'Color',colors(group,:))
            end
        else
            text(xt,yt-offset,['{\it r} = ' num2str(r(2),2) ' (' '{\it p} = ' num2str(p(2),2) ')'],'fontsize',24,'Color',colors(group,:))
        end
    
    end

    %
    ytick = [1:10,20:10:100];
    yticklabel = strings(size(ytick));
    yticklabel(1) = '1'; yticklabel(10) = '10'; yticklabel(19) = '100';
    set(gca,'fontsize',30,...
        'xlim',[-1.2 log10(300)],'xtick',-1:2,'xticklabel',{'0.1','1','10','100'},...
        'ylim',[-10 30],'ytick',-10:10:30,'yticklabel',{'','','','',''})

    switch comparison
        case 1
            text(-1,32,'\bf Sparse Noise','fontsize',24)
            set(gca,'ylim',[-10 30],'ytick',-10:10:30,'yticklabel',{'-10','0','10','20','30'})
            ylabel({'Masking Release','(dB SNR)'})
        case 2
            text(-1,32,'\bf Sparse Noise','fontsize',24)
        case 3
            text(-1,32,'\bf Dense Noise','fontsize',24)
            set(gca,'ylim',[-10 30],'ytick',-10:10:30,'yticklabel',{'-10','0','10','20','30'})
            xlabel({'Frequency Discrimination','(% Difference)'})
            ylabel({'Masking Release','(dB SNR)'})
        case 4
            text(-1,32,'\bf Dense Noise','fontsize',24)
            xlabel({'F0 Discrimination','(% Difference)'})
    end

end

% hack to prevent crashing on mac
drawnow; pause(1);

% save figure
if savefigures
    saveas(f,'figures/fig 8 masking release correlations.png')
end

%% return
return;
%% full correlation table
% This study included six types of measures:
%
% FD : Frequency Discrimination (measured at 1000 Hz)
% F0Ds : Fundamental Frequency Discrimination (measured at 110 Hz)
% SS_F : Stream Segregation with Frequency Differences
% SS_F0 : Stream Segregation with Fundamental Frequency Differences
% SR_Sparse : Speech Reception in Sparse Noise
% SR_Dense : Speech Rception in Dense Noise
%
% In addition, two derived metrics were tabulated:
%
% MR_Sparse : Masking Release in Sparse Noise 
% MR_Dense : Masking Release in Dense Noise
%
% Consequently, there are 8 basic measures. Since the correlation table is
% symmetric and we do not care about self-correlation, there are a total of 
% (N^2 / 2 ) - N relevant correlations ( 8^2 / 2 ) - 8 = 24
%
% Measures were collected for both CI and NH.
% For CI users, the relevant variables are:
% 
% CI_FDT: subject x repetition
% CI_F0DT: subject x repetition
% CI_SS_F: subject x condition x repetition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SS_F0: subject x condition x repetition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SR_Sparse: subject x condition ( condition order: 0, 3, 6, 12, Quiet )
% CI_SR_Dense: subject x condition ( condition order: 0, 3, 6, 12, Quiet )
%
% and likewise for NH.
%
% Make a matrix of summary statistics:
CI_Summary = []; NH_Summary = [];
index = 2:4;
CI_Summary = [
    mean(CI_FDT,2) ...
    mean(CI_F0DT,2) ... 
    mean(mean(CI_SS_F(:,index,:),3),2), ...
    mean(mean(CI_SS_F0(:,index,:),3),2) ...
    mean(mean(CI_SR_Sparse(:,index,:),3),2) ...
    mean(mean(CI_SR_Dense(:,index,:),3),2) ...
    mean(mean(CI_SR_Sparse(:,index,:),3),2) - mean(CI_SR_Sparse(:,1,:),3) ...
    mean(mean(CI_SR_Dense(:,index,:),3),2) - mean(CI_SR_Dense(:,1,:),3) ...
];

NH_Summary = [
    mean(NH_FDT,2) ...
    mean(NH_F0DT,2) ... 
    mean(mean(NH_SS_F(:,index,:),3),2), ...
    mean(mean(NH_SS_F0(:,index,:),3),2) ...
    mean(mean(NH_SR_Sparse(:,index,:),3),2) ...
    mean(mean(NH_SR_Dense(:,index,:),3),2) ...
    mean(mean(NH_SR_Sparse(:,index,:),3),2) - mean(NH_SR_Sparse(:,1,:),3) ...
    mean(mean(NH_SR_Dense(:,index,:),3),2) - mean(NH_SR_Dense(:,1,:),3) ...
];

for a = 1:2
    switch a
        case 1
            X = CI_Summary;
        case 2
            X = NH_Summary;
    end

    disp(' '); disp(a);
    [r,p] = corrcoef(X)
    % m = mean(X); disp(m);
    % s = std(X); disp(s);

end

%% old analyses for listening mode (headphones versus streaming)

% x = CI_Summary(:,1);
% y = CI_Summary(:,2);
% group = categorical(CI_Mode(:,1));
% [h,atab,ctab] = aoctool(x,y,group);
% msgbox(num2str(sqrt(1 - (atab{5,3} / sum([atab{2:5,3}])))))
% corrcoef(x,y)